function [Images] = getFolderPics(FolderPath,Mode,NSequence)
%[Images,Avg] = getFolderPics(FolderPath,Mode,NSequence)
%Modes are 'Single', which returns averages of individual pictures,
%'Normal', which returns averages of OD pictures, and 'Sequential'.
%Sequential mode assumes you've taken repeats of a sequences of NSequence
%pictures, and averages pictures {1, NSequence+1,2*NSequence+1,...} then
%separately {2,NSequence+2,2*NSequence+2,...}
%If an image file contains multiple images, these will be returned as an
%extra dimension in the array. e.g. if there are two images per file, Avg
%will be an array with size [Ny,Nx,Nimg] and Images will be an array with
%size [Ny,Nx,Nimg,Nfiles].
if ~exist('Mode','var')
    Mode = 'Normal';
end


OriginalPath = pwd;
cd(FolderPath)
Files=[dir('*.aia'); dir('*.fits')];

Nfiles=length(Files);
FirstPath = fullfile(FolderPath,Files(1).name);
FirstImgs = readimg(FirstPath);
ImageErrorPositions = [];

if size(FirstImgs,3) == 1 && strcmp(Mode,'Normal')
    %If first image has three pictures, assume all are single OD pics.
    [Vsize,Hsize] = size(FirstImgs);
    Images = zeros(Vsize,Hsize,Nfiles);
    Images(:,:,1) = FirstImgs;
    for ii=2:1:length(Files)
        try
            TmpPath=fullfile(FolderPath,Files(ii).name);
        	TmpImgs = readimg(TmpPath);
            Images(:,:,ii)=TmpImgs;
        catch
            Images(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    [Ny,Nx,Nimg] = size(Images);

elseif size(FirstImgs,3) == 3 && strcmp(Mode,'Normal')
    %If first image has three pictures, assume all are single OD pics.
    FirstOD = getOD(FirstImgs);
    [Vsize,Hsize] = size(FirstOD);
    Images = zeros(Vsize,Hsize,Nfiles);
    Images(:,:,1) = FirstOD;
    for ii=2:1:length(Files)
        try
            TmpPath=fullfile(FolderPath,Files(ii).name);
        	TmpImgs = readimg(TmpPath);
            Images(:,:,ii)=getOD(TmpImgs);
        catch
            Images(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    Avg = sum(Images,3)/size(Images,3);
    [Ny,Nx,Nimg] = size(Images);
    
elseif size(FirstImgs,3) == 4 && strcmp(Mode,'Normal')
    %If first image has three pictures, assume all are single OD pics.
    FirstOD = getOD(FirstImgs);
    [Vsize,Hsize] = size(FirstOD);
    Images = zeros(Vsize,Hsize,Nfiles);
    Images(:,:,1) = FirstOD;
    for ii=2:1:length(Files)
        try
            TmpPath=fullfile(FolderPath,Files(ii).name);
        	TmpImgs = readimg(TmpPath);
            Images(:,:,ii)=getOD(TmpImgs(:,:,[2,3,4]));
        catch
            Images(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    Avg = sum(Images,3)/size(Images,3);
    [Ny,Nx,Nimg] = size(Images);
%     CascadedImage = reshape(Images,[Ny,Nx*Nimg]);
elseif size(FirstImgs,3) == 5 && strcmp(Mode,'Normal')
    %assume double shots. Returns average array of size [Ny,Nx,2] and all
    %images as array of size [Ny,Nx,2,Nimgs]
    FirstOD1 = getOD(FirstImgs(:,:,[2,4,5]));
    FirstOD2 = getOD(FirstImgs(:,:,[3,4,5]));
    [Vsize,Hsize] = size(FirstOD1);
    Images1 = zeros(Vsize,Hsize,Nfiles);
    Images2 = zeros(Vsize,Hsize,Nfiles);
    Images1(:,:,1) = FirstOD1;
    Images2(:,:,1) = FirstOD2;
    for ii=2:1:length(Files)
        try
        TmpPath=fullfile(FolderPath,Files(ii).name);
        TmpImgs = readimg(TmpPath);
        Images1(:,:,ii)=getOD(TmpImgs(:,:,[2,4,5]));
        Images2(:,:,ii) = getOD(TmpImgs(:,:,[3,4,5]));
        catch
            Images1(:,:,ii) = zeros(Vsize,Hsize);
            Images2(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    Images = cat(4,Images1,Images2);
    Avg1 = sum(Images1,3)/size(Images1,3);
    Avg2 = sum(Images2,3)/size(Images2,3);
    Avg = cat(3,Avg1,Avg2);
    [Ny,Nx,Nimg] = size(Images1);
%     CascadedImage1 = reshape(Images1,[Ny,Nx*Nimg]);
%     CascadedImage2 = reshape(Images2,[Ny,Nx*Nimg]);
%     CascadedImage = cat(3,CascadedImage1, CascadedImage2);
    
elseif size(FirstImgs,3) == 6 && strcmp(Mode,'Normal')
    %assume double shots. Returns average array of size [Ny,Nx,2] and all
    %images as array of size [Ny,Nx,2,Nimgs]
    FirstOD1 = getOD(FirstImgs(:,:,[2,4,6]));
    FirstOD2 = getOD(FirstImgs(:,:,[3,5,6]));
    [Vsize,Hsize] = size(FirstOD1);
    Images1 = zeros(Vsize,Hsize,Nfiles);
    Images2 = zeros(Vsize,Hsize,Nfiles);
    Images1(:,:,1) = FirstOD1;
    Images2(:,:,1) = FirstOD2;
    %for ii=4:2:length(Files)
    for ii = 2:1:length(Files)
        try
                TmpPath=fullfile(FolderPath,Files(ii).name);
                TmpImgs = readimg(TmpPath);
                Images1(:,:,ii)=getOD(TmpImgs(:,:,[2,4,6]));
                Images2(:,:,ii) = getOD(TmpImgs(:,:,[3,5,6]));
        catch
            Images1(:,:,ii) = zeros(Vsize,Hsize);
            Images2(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    Images = cat(4,Images1,Images2);
    Avg1 = sum(Images1,3)/size(Images1,3);
    Avg2 = sum(Images2,3)/size(Images2,3);
    Avg = cat(3,Avg1,Avg2);
    [Ny,Nx,Nimg] = size(Images1);
%     CascadedImage1 = reshape(Images1,[Ny,Nx*Nimg]);
%     CascadedImage2 = reshape(Images2,[Ny,Nx*Nimg]);
%     CascadedImage = cat(3,CascadedImage1, CascadedImage2); 
elseif size(FirstImgs,3) == 7 && strcmp(Mode,'Normal')
        %assume double shots. Returns average array of size [Ny,Nx,2] and all
    %images as array of size [Ny,Nx,2,Nimgs]
    FirstOD1 = getOD(FirstImgs(:,:,[2,5,7]));
    FirstOD2 = getOD(FirstImgs(:,:,[3,6,7]));
    [Vsize,Hsize] = size(FirstOD1);
    Images1 = zeros(Vsize,Hsize,Nfiles);
    Images2 = zeros(Vsize,Hsize,Nfiles);
    Images1(:,:,1) = FirstOD1;
    Images2(:,:,1) = FirstOD2;
    %for ii=4:2:length(Files)
    for ii = 2:1:length(Files)
        try
                TmpPath=fullfile(FolderPath,Files(ii).name);
                TmpImgs = readimg(TmpPath);
                Images1(:,:,ii)=getOD(TmpImgs(:,:,[2,5,7]));
                Images2(:,:,ii) = getOD(TmpImgs(:,:,[3,6,7]));
        catch
            Images1(:,:,ii) = zeros(Vsize,Hsize);
            Images2(:,:,ii) = zeros(Vsize,Hsize);
            cat(2,ImageErrorPositions,ii);
        end
    end
    
    Images = cat(4,Images1,Images2);
    Avg1 = sum(Images1,3)/size(Images1,3);
    Avg2 = sum(Images2,3)/size(Images2,3);
    Avg = cat(3,Avg1,Avg2);
    
elseif strcmp(Mode,'Single')%||size(FirstImgs,3)~=3||size(FirstImgs,3)~=5||size(FirstImgs,3)~=6
    %returns average of each individual frame as an array of size
    %[Ny,Nx,Nimg]. Returns all images as array of size [Ny,Nx,Nimg,Nfiles]
    [Vsize,Hsize,Lay] = size(FirstImgs);
    Images = zeros(Vsize,Hsize,Lay,Nfiles);
    Images(:,:,:,1) = FirstImgs;
    %for ii = 4:2:length(Files)
    for ii = 2:1:length(Files)
        try
            TmpPath=fullfile(FolderPath,Files(ii).name);
            TmpImgs = readimg(TmpPath);
            Images(:,:,:,ii) = TmpImgs;
        catch
            Images(:,:,:,ii) = zeros(Vsize,Hsize,Lay);
            cat(2,ImageErrorPositions,ii);
        end
    end
    Avg = sum(Images,4)/size(Images,4);
    %CascadedImage = 0;
    
elseif strcmp(Mode,'Sequential')
    [Vsize,Hsize,Lay] = size(FirstImgs);
    Images = zeros(Vsize,Hsize,Lay,Nfiles);
    Images(:,:,:,1) = FirstImgs;
    CurrentSequence = 1;
    for ii = 2:1:length(Files)
        CurrentSequence = mod(CurrentSequence+1,NSequence);
        try
            TmpPath=fullfile(FolderPath,Files(ii).name);
            TmpImgs = readimg(TmpPath);
            Images(:,:,:,ii) = TmpImgs;
        catch
            Images(:,:,:,ii) = zeros(Vsize,Hsize,Lay);
            cat(2,ImageErrorPositions,ii);
        end
    end
    NPartialSequences = ceil(length(Files)/NSequence);
    SequenceCounter = kron(ones(1,NPartialSequences),1:NSequence);
    SequenceCounter = SequenceCounter(1:length(Files));
    Avg = zeros(Vsize,Hsize,Lay,NSequence);
    for jj = 1:NSequence
        CurrentSeqImages = Images(:,:,:,SequenceCounter==jj);
        Avg(:,:,:,jj) = sum(CurrentSeqImages,4)/size(CurrentSeqImages,4);
    end
        
        
end
cd(OriginalPath);

end

